import osmnx as ox
import networkx as nx
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
Create a graph from a place¶
# Using osmnx.graph_from_address to turn textual description of place to networkx graph of that place map
place_name = "university of toronto"
graph = ox.graph_from_address(place_name)
# graph = ox.graph_from_place(place_name, network_type='drive') # retrieve drive netowrk only
Visualizing the graph¶
# Visualize that graph using matplotlib
ox.plot_graph(graph, figsize=(10,10))
(<Figure size 720x720 with 1 Axes>, <AxesSubplot:>)
# Let's prune the graph a bit by changing the width of our scope
graph = ox.graph_from_address("university of toronto", dist = 300)
# Visualize graph nodes with folium
ox.folium.plot_graph_folium(graph)
# add markers for points of interest
import folium
m = folium.Map(location=[43.663774, -79.395981], tiles="OpenStreetMap", zoom_start=16)
center=(43.662643, -79.395689) # UofT main building
source_point = (43.664527, -79.392442) # King Edward VII Equestrian Statue
destination_point = (43.659659, -79.397669) # Bahen Centre for Information Technology at UofT
folium.Marker(location = [center[0], center[1]], popup = 'UofT main building', icon=folium.Icon(color='blue', icon='star', prefix='fa')).add_to(m)
folium.Marker(location = [source_point[0],source_point[1]], popup = 'King Edward VII Equestrian Statue', icon = folium.Icon(color='red', icon='home', prefix='fa')
).add_to(m)
folium.Marker(location = [destination_point[0],destination_point[1]], popup = 'Bahen Centre for Information Technology at UofT', icon = folium.Icon(color='green', icon='flag', prefix='fa')
).add_to(m)
m
center[0]
43.662643
# plotting map using ipyleaflet
from ipyleaflet import *
center=(43.662643, -79.395689) # UofT main building
source_point = (43.664527, -79.392442) # King Edward VII Equestrian Statue
destination_point = (43.659659, -79.397669) # Bahen Centre for Information Technology at UofT
m = Map(center=center, zoom=15)
marker = Marker(location=center, draggable=False)
# King Edward VII Equestrian Statue
marker1 = Marker(name='marker1', location=source_point, icon=Icon(icon_url='https://raw.githubusercontent.com/pointhi/leaflet-color-markers/master/img/marker-icon-red.png', icon_size=[25, 41], icon_anchor=[12, 41]))
# Bahen Centre for Information Technology
marker2 = Marker(name='marker2', location=destination_point, icon=Icon(icon_url='https://raw.githubusercontent.com/pointhi/leaflet-color-markers/master/img/marker-icon-yellow.png', icon_size=[25, 41], icon_anchor=[12, 41]))
# for more icons: https://ipyleaflet.readthedocs.io/en/latest/api_reference/awesome_icon.html
m.add_layer(marker1)
m.add_layer(marker2)
control = LayersControl(position='topright')
m.add_control(control)
m.add_layer(marker);
m
Extracting graph information¶
# type of the generated graph
type(graph)
networkx.classes.multidigraph.MultiDiGraph
# extract the nodes and edges from the graph as GeoDataFrames
nodes, edges = ox.graph_to_gdfs(graph)
nodes.head(5)
| y | x | osmid | highway | geometry | |
|---|---|---|---|---|---|
| 6391301132 | 43.662454 | -79.399407 | 6391301132 | NaN | POINT (-79.39941 43.66245) |
| 1258707987 | 43.662295 | -79.395994 | 1258707987 | NaN | POINT (-79.39599 43.66229) |
| 1258707990 | 43.661657 | -79.396219 | 1258707990 | NaN | POINT (-79.39622 43.66166) |
| 389678111 | 43.662654 | -79.394529 | 389678111 | NaN | POINT (-79.39453 43.66265) |
| 389678112 | 43.662613 | -79.394646 | 389678112 | NaN | POINT (-79.39465 43.66261) |
edges.head(5)
| osmid | highway | access | oneway | length | geometry | lanes | name | maxspeed | service | tunnel | u | v | key | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 234732970 | footway | permissive | False | 10.945 | LINESTRING (-79.39941 43.66245, -79.39954 43.6... | NaN | NaN | NaN | NaN | NaN | 6391301132 | 2428750571 | 0 |
| 1 | 682466783 | footway | NaN | False | 7.410 | LINESTRING (-79.39941 43.66245, -79.39943 43.6... | NaN | NaN | NaN | NaN | NaN | 6391301132 | 6391312709 | 0 |
| 2 | [234732970, 396822863, 366318543, 639850036, 2... | [corridor, footway, steps] | permissive | False | 111.947 | LINESTRING (-79.39941 43.66245, -79.39935 43.6... | NaN | NaN | NaN | NaN | NaN | 6391301132 | 6028561930 | 0 |
| 3 | 34011873 | footway | NaN | False | 14.600 | LINESTRING (-79.39599 43.66229, -79.39582 43.6... | NaN | NaN | NaN | NaN | NaN | 1258707987 | 7311083158 | 0 |
| 4 | 204343740 | footway | permissive | False | 13.855 | LINESTRING (-79.39599 43.66229, -79.39604 43.6... | NaN | NaN | NaN | NaN | NaN | 1258707987 | 2143485199 | 0 |
# Retrieve only edges from the graph
edges = ox.graph_to_gdfs(graph, nodes=False, edges=True)
# rendering the 2nd node
list(graph.nodes(data=True))[1]
(1258707987, {'y': 43.6622948, 'x': -79.3959935, 'osmid': 1258707987})
# rendering the 1std edge
list(graph.edges(data=True))[1]
(6391301132,
6391312709,
{'osmid': 682466783, 'highway': 'footway', 'oneway': False, 'length': 7.41})
# get road types
print(edges['highway'].value_counts())
footway 836
service 106
residential 62
path 60
tertiary 54
[footway, steps] 18
unclassified 17
pedestrian 10
steps 4
[service, footway] 4
[path, steps] 4
[corridor, footway, steps] 2
secondary 2
[footway, path] 2
[service, path, steps] 2
Name: highway, dtype: int64
# Converting node and edge GeoDataFrames to a MultiDiGraph.
graph_2=ox.graph_from_gdfs(nodes, edges)
ig2, ax2 = ox.plot_graph(graph_2,figsize=(10,10))
# Calculate network statistics
stats = ox.basic_stats(graph, circuity_dist='euclidean')
stats
{'n': 422,
'm': 1183,
'k_avg': 5.606635071090047,
'intersection_count': 366,
'streets_per_node_avg': 3.059241706161137,
'streets_per_node_counts': {0: 0, 1: 56, 2: 0, 3: 234, 4: 128, 5: 3, 6: 1},
'streets_per_node_proportion': {0: 0.0,
1: 0.13270142180094788,
2: 0.0,
3: 0.5545023696682464,
4: 0.3033175355450237,
5: 0.0071090047393364926,
6: 0.002369668246445498},
'edge_length_total': 33526.35599999993,
'edge_length_avg': 28.340114961961053,
'street_length_total': 17284.822999999993,
'street_length_avg': 28.664714759535645,
'street_segments_count': 603,
'node_density_km': None,
'intersection_density_km': None,
'edge_density_km': None,
'street_density_km': None,
'circuity_avg': 95994.84522129643,
'self_loop_proportion': 0.0,
'clean_intersection_count': None,
'clean_intersection_density_km': None}
# Get the Convex Hull of the network
convex_hull = edges.unary_union.convex_hull
# Show output
convex_hull
#extended statistics for the network using the Convex Hull
# Calculate the area
area = convex_hull.area
# Calculate statistics with density information
stats = ox.basic_stats(graph, area=area)
extended_stats = ox.extended_stats(graph, ecc=True, cc=True)
# Add extened statistics to the basic statistics
for key, value in extended_stats.items():
stats[key] = value
# Convert the dictionary to a Pandas series for a nicer output
pd.Series(stats)
n 422
m 1183
k_avg 5.60664
intersection_count 366
streets_per_node_avg 3.05924
streets_per_node_counts {0: 0, 1: 56, 2: 0, 3: 234, 4: 128, 5: 3, 6: 1}
streets_per_node_proportion {0: 0.0, 1: 0.13270142180094788, 2: 0.0, 3: 0....
edge_length_total 33526.4
edge_length_avg 28.3401
street_length_total 17284.8
street_length_avg 28.6647
street_segments_count 603
node_density_km 1.1151e+13
intersection_density_km 9.67124e+12
edge_density_km 8.85906e+14
street_density_km 4.56737e+14
circuity_avg 1.0287
self_loop_proportion 0
clean_intersection_count None
clean_intersection_density_km None
avg_neighbor_degree {6391301132: 3.3333333333333335, 1258707987: 2...
avg_neighbor_degree_avg 3.05458
avg_weighted_neighbor_degree {6391301132: 0.07674479286580405, 1258707987: ...
avg_weighted_neighbor_degree_avg 0.179252
degree_centrality {6391301132: 0.014251781472684086, 1258707987:...
degree_centrality_avg 0.0133174
clustering_coefficient {6391301132: 0, 1258707987: 0, 1258707990: 0, ...
clustering_coefficient_avg 0.0511848
clustering_coefficient_weighted {6391301132: 0, 1258707987: 0, 1258707990: 0, ...
clustering_coefficient_weighted_avg 0.00866079
pagerank {6391301132: 0.00272874173708308, 1258707987: ...
pagerank_max_node 2143436415
pagerank_max 0.00520859
pagerank_min_node 8544407502
pagerank_min 0.000394535
eccentricity {6391301132: 716.754, 1258707987: 745.887, 125...
diameter 1119.22
radius 555.92
center [389677904]
periphery [24960076]
closeness_centrality {6391301132: 0.002814157839676847, 1258707987:...
closeness_centrality_avg 0.00252743
dtype: object
Projection¶
# Check crs
edges.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
# # Reproject to Mercator
# edges_proj = edges_proj.to_crs(epsg=3857)
# edges_proj.crs
Shortest path analysis¶
# getting the nearest nodes to the locations
orig_node = ox.get_nearest_node(graph, source_point)
target_node = ox.get_nearest_node(graph, destination_point)
# Render the closest node id to origin and destination points
target_node, target_node
(239055725, 239055725)
# Retrieve the rows from the nodes GeoDataFrame
o_closest = nodes.loc[orig_node]
t_closest = nodes.loc[target_node]
# Create a GeoDataFrame from the origin and target points
od_nodes = gpd.GeoDataFrame([o_closest, t_closest], geometry='geometry', crs=nodes.crs)
od_nodes.head()
| y | x | osmid | highway | geometry | |
|---|---|---|---|---|---|
| 390545921 | 43.665132 | -79.394109 | 390545921 | NaN | POINT (-79.39411 43.66513) |
| 239055725 | 43.660783 | -79.397878 | 239055725 | NaN | POINT (-79.39788 43.66078) |
# Calculate the shortest path
route = nx.shortest_path(G=graph, source=orig_node, target=target_node, weight='length')
# Show what we have
print(route)
[390545921, 60654129, 60654119, 389678001, 389678002, 2143434369, 390550470, 127289393, 8277128565, 8277128566, 4920594801, 3996671922, 80927418, 7967019552, 127275360, 80927426, 2143468197, 2143468182, 55808564, 55808527, 130170945, 389677905, 389678182, 389677906, 50885141, 389678180, 1258706668, 2143436407, 1258706673, 2143402269, 2143402268, 1258706670, 239055725]
# Plot the shortest path
fig, ax = ox.plot_graph_route(graph, route, figsize=(15,15))
# """Draw leaflet map based for the graph passed to the function
# with highlighting certain nodes with a marker. ipyleaflet rendering
# can be a little slow when dealing with graphs with many node, in that
# case when number of nodes > 1000, the function reverts to folium for rendering
# which is a plain black and white rendering of the map not like ipyleaflet
# Parameters
# ----------
# G: networkx graph returned by osmnx
# highlight: list of nodes to be marked on the leaflet map
# zoom: the initial zooming level for the rendered map
# Returns
# -------
# m: ipyleaflet/folium map for the osmnx graph with optional highlighted nodes
# """
import ipyleaflet as lf
def draw_route(G, route, zoom = 16):
if len(G) >= 1000:
print(f"The graph has {len(G)} which is a lot, we will use basic faster folium instead")
m = ox.plot_route_folium(G = G, route = route)
return m
center_osmid = ox.stats.extended_stats(G,ecc=True)['center'][0]
G_gdfs = ox.graph_to_gdfs(G)
nodes_frame = G_gdfs[0]
ways_frame = G_gdfs[1]
center_node = nodes_frame.loc[center_osmid]
location = (center_node['y'], center_node['x'])
m = lf.Map(center = location, zoom = zoom)
start_node = nodes_frame.loc[route[0]]
end_node = nodes_frame.loc[route[len(route)-1]]
start_xy = (start_node['y'], start_node['x'])
end_xy = (end_node['y'], end_node['x'])
marker = lf.Marker(location = start_xy, draggable = False)
m.add_layer(marker)
marker = lf.Marker(location = end_xy, draggable = False)
m.add_layer(marker)
for u, v in zip(route[0:], route[1:]):
try:
x, y = (ways_frame.query(f'u == {u} and v == {v}').to_dict('list')['geometry'])[0].coords.xy
except:
x, y = (ways_frame.query(f'u == {v} and v == {u}').to_dict('list')['geometry'])[0].coords.xy
points = map(list, [*zip([*y],[*x])])
ant_path = lf.AntPath(
locations = [*points],
dash_array=[1, 10],
delay=1000,
color='red',
pulse_color='black'
)
m.add_layer(ant_path)
return m
draw_route(graph, route)
Retrieve buildings information¶
buildings = ox.footprints_from_place(place_name)
buildings_proj = buildings.to_crs(crs=edges.crs)
buildings_proj.head(5)
| nodes | addr:housenumber | addr:street | building | building:colour | building:material | college | name | operator | geometry | ... | name:ne | shop | source | roof:levels | wikimedia_commons | building:min_level | members | opening_hours | type | addr:full | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7685383 | [55808755, 340165033, 55808760, 55808764, 8504... | 300 | Huron Street | residential | orange | brick | New College | Wetmore Hall | University of Toronto | POLYGON ((-79.39992 43.66229, -79.40009 43.662... | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 7948423 | [59227761, 2143460128, 3702927984, 2143460131,... | 100 | St George Street | university | NaN | NaN | NaN | Sidney Smith Hall | University of Toronto | POLYGON ((-79.39948 43.66268, -79.39942 43.662... | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 7991747 | [59817417, 59817411, 59817412, 3033737068, 303... | 130 | St George Street | university | gray | concrete | NaN | Robarts Library | University of Toronto | POLYGON ((-79.39972 43.66396, -79.39972 43.664... | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 8027952 | [60006599, 1633499537, 1633499539, 60006598, 8... | 113 | St. Joseph Street | university | gray | concrete | St. Michael's College | J. M. Kelly Library | University of Toronto | POLYGON ((-79.38912 43.66470, -79.38917 43.664... | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 8027962 | [60006651, 2143468180, 60006652, 60006653, 600... | 73-75 | St George Street | university | NaN | NaN | University College | Sir Daniel Wilson Residence (University College) | University of Toronto | POLYGON ((-79.39747 43.66297, -79.39755 43.662... | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
5 rows × 57 columns
# Get the nodes along the shortest path
route_nodes = nodes.loc[route]
route_nodes
| y | x | osmid | highway | geometry | |
|---|---|---|---|---|---|
| 390545921 | 43.665132 | -79.394109 | 390545921 | NaN | POINT (-79.39411 43.66513) |
| 60654129 | 43.665060 | -79.394333 | 60654129 | NaN | POINT (-79.39433 43.66506) |
| 60654119 | 43.664991 | -79.394661 | 60654119 | NaN | POINT (-79.39466 43.66499) |
| 389678001 | 43.664946 | -79.394875 | 389678001 | NaN | POINT (-79.39487 43.66495) |
| 389678002 | 43.664826 | -79.395444 | 389678002 | NaN | POINT (-79.39544 43.66483) |
| 2143434369 | 43.664756 | -79.395777 | 2143434369 | NaN | POINT (-79.39578 43.66476) |
| 390550470 | 43.664749 | -79.395809 | 390550470 | NaN | POINT (-79.39581 43.66475) |
| 127289393 | 43.664672 | -79.395777 | 127289393 | traffic_signals | POINT (-79.39578 43.66467) |
| 8277128565 | 43.664601 | -79.395749 | 8277128565 | NaN | POINT (-79.39575 43.66460) |
| 8277128566 | 43.664597 | -79.395769 | 8277128566 | NaN | POINT (-79.39577 43.66460) |
| 4920594801 | 43.664465 | -79.395717 | 4920594801 | NaN | POINT (-79.39572 43.66447) |
| 3996671922 | 43.663604 | -79.395370 | 3996671922 | NaN | POINT (-79.39537 43.66360) |
| 80927418 | 43.663569 | -79.395356 | 80927418 | NaN | POINT (-79.39536 43.66357) |
| 7967019552 | 43.663420 | -79.395484 | 7967019552 | NaN | POINT (-79.39548 43.66342) |
| 127275360 | 43.663258 | -79.396201 | 127275360 | NaN | POINT (-79.39620 43.66326) |
| 80927426 | 43.663136 | -79.396929 | 80927426 | NaN | POINT (-79.39693 43.66314) |
| 2143468197 | 43.663041 | -79.396876 | 2143468197 | NaN | POINT (-79.39688 43.66304) |
| 2143468182 | 43.662965 | -79.396830 | 2143468182 | NaN | POINT (-79.39683 43.66296) |
| 55808564 | 43.662909 | -79.396797 | 55808564 | NaN | POINT (-79.39680 43.66291) |
| 55808527 | 43.662641 | -79.396999 | 55808527 | NaN | POINT (-79.39700 43.66264) |
| 130170945 | 43.662550 | -79.397394 | 130170945 | NaN | POINT (-79.39739 43.66255) |
| 389677905 | 43.662495 | -79.397650 | 389677905 | NaN | POINT (-79.39765 43.66250) |
| 389678182 | 43.662183 | -79.397533 | 389678182 | NaN | POINT (-79.39753 43.66218) |
| 389677906 | 43.662069 | -79.397487 | 389677906 | NaN | POINT (-79.39749 43.66207) |
| 50885141 | 43.662049 | -79.397593 | 50885141 | crossing | POINT (-79.39759 43.66205) |
| 389678180 | 43.662033 | -79.397663 | 389678180 | NaN | POINT (-79.39766 43.66203) |
| 1258706668 | 43.661075 | -79.397266 | 1258706668 | NaN | POINT (-79.39727 43.66107) |
| 2143436407 | 43.661055 | -79.397365 | 2143436407 | NaN | POINT (-79.39737 43.66105) |
| 1258706673 | 43.660930 | -79.397436 | 1258706673 | NaN | POINT (-79.39744 43.66093) |
| 2143402269 | 43.660887 | -79.397637 | 2143402269 | NaN | POINT (-79.39764 43.66089) |
| 2143402268 | 43.660862 | -79.397754 | 2143402268 | NaN | POINT (-79.39775 43.66086) |
| 1258706670 | 43.660832 | -79.397896 | 1258706670 | NaN | POINT (-79.39790 43.66083) |
| 239055725 | 43.660783 | -79.397878 | 239055725 | NaN | POINT (-79.39788 43.66078) |
from shapely.geometry import LineString, Point
# Create a geometry for the shortest path
route_line = LineString(list(route_nodes.geometry.values))
route_line
# Create a GeoDataFrame
route_geom = gpd.GeoDataFrame([[route_line]], geometry='geometry', crs=edges.crs, columns=['geometry'])
# Add a list of osmids associated with the route
route_geom.loc[0, 'osmids'] = str(list(route_nodes['osmid'].values))
# Calculate the route length
route_geom['length_m'] = route_geom.length
route_geom.head()
C:\Anaconda3\envs\uoft\lib\site-packages\ipykernel_launcher.py:8: UserWarning: Geometry is in a geographic CRS. Results from 'length' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
| geometry | osmids | length_m | |
|---|---|---|---|
| 0 | LINESTRING (-79.39411 43.66513, -79.39433 43.6... | [390545921, 60654129, 60654119, 389678001, 389... | 0.00841 |
# Plot edges and nodes
ax = edges.plot(linewidth=0.75, color='gray', figsize=(15,15))
ax = nodes.plot(ax=ax, markersize=2, color='gray')
# Add buildings
ax = buildings_proj.plot(ax=ax, facecolor='khaki', alpha=0.7)
# Add the route
ax = route_geom.plot(ax=ax, linewidth=2, linestyle='--', color='red')
# Add the origin and destination nodes of the route
ax = od_nodes.plot(ax=ax, markersize=24, color='green')
Finding points of interest (POIs) using overpass-api¶
import requests
import json
overpass_url = "http://overpass-api.de/api/interpreter"
overpass_query = """
[out:json];
(node["amenity"="restaurant"](around:1000,43.66, -79.39);
way["amenity"="restaurant"](around:1000,43.66, -79.39);
rel["amenity"="restaurant"](around:1000,43.66, -79.39);
);
out center;
"""
response = requests.get(overpass_url, params={'data': overpass_query})
restaurant = response.json()
len(restaurant['elements'])
196
# Extract the lon, lat and name of each restaurant:
lons = []
lats = []
text = []
for elem in restaurant['elements']:
if elem['type'] == 'node':
lons.append(elem['lon'])
lats.append(elem['lat'])
elif 'center' in elem:
lons.append(elem['center']['lon'])
lats.append(elem['center']['lat'])
else: pass
if 'name' not in elem['tags']:
text.append('Noname')
else:
text.append(elem['tags']['name'])
restaurants = dict(type='scattermapbox',
lat=lats,
lon=lons,
mode='markers',
text=text,
marker=dict(size=8, color='blue'),
hoverinfo='text',
showlegend=False)
# plotting restaurants' locations around University of Toronto
import plotly.graph_objects as obj
center=(43.662643, -79.395689) # UofT main building
fig = obj.Figure(obj.Scattermapbox(restaurants))
# defining plot layout
fig.update_layout(mapbox_style="stamen-terrain")
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0}, mapbox = {'center': {'lat': center[0], 'lon': center[1]}, 'zoom': 13})
fig.show()